Identification of PIK3CG as a hub in septic myocardial injury using network pharmacology and weighted gene co‐expression network analysis

Abstract Sepsis causes multiple organ injuries, among which the heart is one most severely damaged organ. Melatonin (MEL) alleviates septic myocardial injury, although a systematic and comprehensive approach is still lacking to understand the precise protective machinery of MEL. This study aimed to examine the underlying mechanisms of MEL on improvement of septic myocardial injury at a systematic level. This study integrated three analytic modalities including database investigations, RNA‐seq analysis, and weighted gene co‐expression network analysis (WCGNA), in order to acquire a set of genes associated with the pathogenesis of sepsis. The Drugbank database was employed to predict genes that may serve as pharmacological targets for MEL‐elicited benefits, if any. A pharmacological protein–protein interaction network was subsequently constructed, and 66 hub genes were captured which were enriched in a variety of immune response pathways. Notably, PIK3CG, one of the hub genes, displayed high topological characteristic values, strongly suggesting its promise as a novel target for MEL‐evoked treatment of septic myocardial injury. Importantly, molecular docking simulation experiments as well as in vitro and in vivo studies supported an essential role for PIK3CG in MEL‐elicited effect on septic myocardial injury. This study systematically clarified the mechanisms of MEL intervention in septic myocardial injury involved multiple targets and multiple pathways. Moreover, PIK3CG‐governed signaling cascade plays an important role in the etiology of sepsis and septic myocardial injury. Findings from our study provide valuable information on novel intervention targets for the management of septic myocardial injury. More importantly, this study has indicated the utility of combining a series of techniques for disease target discovery and exploration of possible drug targets, which should shed some light on elucidation of experimental and clinical drug action mechanisms systematically.

. Sepsis is believed the main culprit for high morbidity and mortality in the intensive care unit. 1 In spite of recent advances in intensive care and the discovery of antibiotics, the mortality rate remains high for severe sepsis, with an everrising incidence rate each year. 2 Among various organs afflicted by sepsis, the heart is deemed a primary target for sepsis-induced MODS. Accumulating evidence has denoted the culprit role for myocardial dysfunction in the deterioration of sepsis complications. 3 The mortality of patients with concurrent sepsis and myocardial dysfunction is approximately 70%-90%, much higher than those patients without myocardial injury (20%). 4 Therefore, it is pertinent to identify novel and effective therapeutic strategies for septic myocardial injury to improve the prognosis of sepsis.
Melatonin (MEL) is a neuroendocrine hormone synthesized by the pineal gland, with an inherent role in circadian rhythm control. 5 In addition, MEL is produced in most organs and tissues at much higher levels than the pineal gland-the so-called extrapineal MEL. 6 The high levels of extrapineal MEL are closely associated with the ability to eliminate reactive oxygen species (ROS), evoke antioxidant enzyme expression, reduce peroxidase activity, and maintain mitochondrial homeostasis, thereby exerting a powerful anti-inflammatory and antioxidant effect. 7,8 Multiple studies have denoted an active role for MEL in various cardiovascular diseases such as hypertension, myocardial ischemia-reperfusion injury, and atherosclerosis. 9 Notably, MEL also interacts with various pleiotropic pathways, including NOD-like receptor family pyrin domain containing 3 (NLRP3) signaling, 10,11 and the mammalian STE20-like kinase-1 (Mst1)/c-Jun N-terminal kinase (JNK) pathway, 12 and was suggested to protect against septic myocardial injury.
Although several studies have shown a protective role for MEL in sepsis and septic myocardial injury, the current understanding with regards to pharmacological actions behind MEL is mainly derived from its intervention with common molecules or signaling pathways. A more thorough, systematic and comprehensive approach is urgently needed to unveil the full spectrum behind MEL-elicited protective mechanisms. Network pharmacology provides a practical strategy for elucidation of potential mechanisms of multicomponent and multitarget drugs by analyzing multiple networks of complex interactions, commonly performed to discover the key targets for disease treatment. 13,14 Moreover, the weighted gene co-expression network analysis (WGCNA) and RNA-seq are also widely used to identify core targets of human diseases. 15,16 Hence, three strategies (RNA-seq analysis, WGCNA, 17 and database search) were conducted in this study to obtain genes associated with the pathogenesis of sepsis ( Figure 1a). Then, pharmacological targets for MEL were predicted according to the Drugbank database. Subsequently, a pharmacological network 13,14 with genes associated with sepsis and putative targets for MEL was constructed to dissect possible therapeutic hub genes.    Figure S1) and raw expression data production. Principal component analysis (PCA) results based on these genes showed that the 12 biological replicates from each group were closely clustered together, indicating an acceptable within group variation in replicates ( Figure 1b). Moreover, the PCA score plot displayed an observable distinction between CLP and Sham groups, revealing that sepsis caused significant changes in gene expression ( Figure 1b). As illustrated in Figure 1c, a total of 778 differentially expressed genes (DEGs) were detected using Deseq2, of which Saa3, 18 Csf3, 19 and 2.2 | Identification of core genes related to SV, CO, sepsis score, and group in the weighted co-expression network WGCNA is generally used to identify modules and candidate marker genes related to external sample traits and has been successfully performed in various biological entities, such as cancer and mouse genome. 21,22 Here in this study, WGCNA was applied to quest genes related with biometric, plasma and echocardiographic parameters of septic mice (see Section 5 for details), in an effort to identify core candidate targets by structuring a scale-free gene co-expression network.
The normalized count matrix of a total of 5699 genes after quality control was used as input data. The power of β = 8 (scale free   Table S2). This study predicted a total of 109 drug targets ( Figure 3a and Table S3) (Table S5). These aforementioned 66 hub genes were enriched into three functional groups (BP, CC, and MF) of GO terms ( Figure 3f and Table S6; only displayed the top 10 terms in each group). Hub genes in the BP group were mainly enriched in response to lipopolysaccharide, response to molecule of bacterial origin, and positive regulation of cytokine production ( Figure 3f, Table S6). The genes in the CC F I G U R E 2 The core genes of WGCNA modules associated with sepsis. (a) The left figure shows the network topology analysis of different soft-thresholding power β. The weighted coefficient β is finally set to 8 (R 2 = 0.85) to construct a signed scale-free network. The figure on the right illustrates the mean network connectivity under different weighted coefficient. (b) Heatmap of the correlation coefficients between modules and biometric, plasma and echocardiographic parameters of septic mice. Each column refers to a parameter, and each row refers to a ME (the first principal component of each gene module). According to the legend on the right, the matrix is color-coded by correlation (red, positive correlation; blue, negative correlation). The Pearson correlation coefficient is shown in the rectangle and the P value is shown in parentheses. Module turquoise is significantly correlated with parameters SV, CO, group and sepsis score. (c) According to the threshold (jMMj > 0.8, jGSj > 0.8; details in Section 5), the core genes correlated with SV, CO, group, and sepsis score are extracted from the module turquoise. The veen diagram of these core genes is drawn in (d), and the GO terms enrichment result is shown in (e). AST, aspartate aminotransferase; BUN, blood urea nitrogen; CO, cardiac output; CK, creatine kinase; GRA, granulocytes; GS, Gene significance; LDH, lactic dehydrogenase; LYM, lymphocytes; MEs, module eigengenes; MID, middle cells; MM, module membership; PLT, platelets; SV, stroke volume; RBC, red blood cells; WBC, white blood cells; WGCNA, weighted gene co-expression network analysis F I G U R E 3 The molecular basis of the protective effect of MEL on sepsis. (a) Venn diagram depicting the overlap between sepsisrelated targets and MEL putative targets. The background image shows the chemical structure of MEL. (b) The molecular docking simulation of the binding pattern of MEL with ADA, PIK3CG, and MAPKAPK3. MEL, melatonin; IL-6, Interleukin 6. (c) The interaction network describes the connections between MEL putative targets (nodes in the left circle) and sepsis-related genes (nodes in the right concentric circle). The nodes circled by dotted lines are genes shared by sepsis-related genes and MEL putative targets (MAPKAPK3 has no link). Yellow and blue nodes symbolize hub genes and non-hub genes, respectively. Node size is ranked according to its network centrality, and PIK3CG is the second largest node in the network, behind the IL-6. The red triangle nodes refer to the top five KEGG pathways (sorted by rich factor) of the signal transduction classification pathways enriched by the hub genes. The bar graph (d) shows the classification of KEGG pathways enriched by these hub genes, the dot plot (e) shows the pathways of signal transduction classification, and the bar chart (f) exhibits the enriched GO terms group were mainly focused on membrane raft and secretory granule lumen, while the genes in the MF group were heavily resided in cytokine activity, immune receptor activity, receptor ligand activity, and signaling receptor activator activity ( Figure 3f, Table S6).
Moreover, the enriched KEGG pathways were classified, which demonstrated an essential involvement in signal transduction and immune response (Figure 3d and Table S7), including PI3K-AKT, Rap1 signaling, NF-kappa B signaling, and FoxO signaling pathways Table S8). Functional analysis of these hub genes suggested that MEL ameliorated myocardial injury in sepsis by multiple targets and pathways.
Among these hub genes, PIK3CG with high centrality values, participated in PI3K-AKT signaling pathway, and cGMP-PKG signaling pathway (Figure 3c, Tables S5 and S8). In addition, PIK3CG showed the strongest binding affinity with MEL among the three common targets shared between MEL targets and sepsis targets ( Figure 3b).
Hence, we propose that MEL plays a protective role in sepsis by regulating immune-related pathways, mainly through the PIK3CG-based PI3K-AKT pathway.

|
The effects of MEL on survival rate, sepsis score, and anal temperature in CLP-injured mice showed that CLP led to a high mortality rate (nearly 75%). However, MEL overtly improved the survival rate. Among various MEL groups, 30 mg/kg MEL pretreatment exhibited an obvious protective effect (vs. CLP group, Figure 4d; p < 0.05). Therefore, this level of MEL was used for the remaining studies. Subsequently, sepsis score and anal temperature were monitored 8 h post-CLP. In this study, the MSS system was introduced to monitor experimental mice based on their physical appearance, level of consciousness, activity, response to stimuli, eyes, respiration rate and quality (from 0 to 4 points for each criteria). 25 The results showed that CLP significantly increased sepsis score, the effects of which were reversed by MEL treatment (Figure 4e). As shown in Figure 4f, MEL significantly increased the anal temperature in CLP-injured mice (vs. CLP group, p < 0.05).
2.5 | The effects of MEL on blood routine parameters, blood biochemical parameters, cardiac function, and myocardial structure in CLP-injured mice

| DISCUSSION
In this study, we employed a combination of RNA-seq analysis, WGCNA, pharmacological network construction, and molecular docking to identify genes and pathways involved in the treatment of septic myocardial injury. First, the 453 sepsis-related genes were obtained by merging the top 200 DEGs from RNA-seq analysis, 109 known targets from database, and 156 core genes from WGCNA. Second, 109 targets of MEL were predicted using the Drugbank database.
Based on a pharmacological network constructed by drug targets and disease genes, PIK3CG was identified as one essential candidate target as supported by molecular docking. As expected, in vivo and in vitro studies also confirmed that MEL exerts a pronounced cardioprotective effect against septic myocardial injury through inhibition of PIK3CG-related signaling pathway (Figure 7).
PI3K is an important component of the PI3K-Akt intracellular signaling pathway, consisting of regulatory subunit (p85) and catalytic subunit (p110). PI3K is reported to be involved in cell proliferation, apoptosis, migration, and inflammation among other pathophysiological processes, affecting the onset and development of a variety of cardiovascular diseases. 26 There are four types of catalytic subunits of PI3K, namely p110α, β, δ, and γ. Among them, p110α, β, and δ are involved in immunity and inflammation. 27,28 PI3Kγ/p110γ, also known as PIK3CG, is a key regulator molecule in the pathological process of inflammation and oxidation, and thus plays a critical role in a variety of cardiovascular diseases. 26 For example, blood pressure can be declined in a dose-dependent fashion through inhibition of PI3KCG in mice. 29 Moreover, the deficiency of PIK3CG protects anthracyclineinduced cardiotoxicity and decreases tumor growth. 30 However, whether PIK3CG participates in septic myocardial injury is essentially unknown. Hence, we performed the network pharmacological analysis and molecular docking simulation and found that PIK3CG was a potential core target of MEL for the treatment of septic myocardial injury. Importantly, in vivo and in vitro experiments showed that MEL alleviated sepsis and septic myocardial injury by inhibiting PIK3CGrelated signaling. Hence, we present evidence here for optimization of a series of methods for discovering disease and drug targets, to identify and verify involvement of PIK3CG in MEL-offered benefit against sepsis myocardial injury, thus unveiling a key role for PI3K family members in sepsis.
As mentioned above, a systematic and comprehensive approach is urgently needed to explore the mechanisms of MEL in the protection against septic myocardial injury. Network pharmacology provides a practical strategy for elucidating the potential mechanisms of multicomponent and multitarget agents by analyzing various networks of complex interactions, which is widely used to discover key pathways and target molecules for disease treatment. 13 However, target information obtained by database retrieval is limited, which compromises the accuracy and comprehensiveness of the mechanisms research. 13 F I G U R E 5 The effects of MEL on myocardial structure, inflammatory response, oxidative stress, and PIK3CG-related signaling pathway in CLP-injured mice. (l) qRT-PCR analysis of PIK3CG, Pde4a, Syk, and Myc mRNA by normalizing to β-actin. *p < 0.05, **p < 0.01, and ***p < 0.001 compared with the CLP group; ns, nonsignificant; n = 6. Using one-way ANOVA and Tukey's multiple comparisons. CLP, cecal ligation and puncture; H&E, hematoxylin-eosin; IHC, immunohistochemistry; IL-6, interleukin 6; MEL, melatonin; NLRP3, NOD-like receptor family pyrin domain containing 3; TNF-a, tumor necrosis factor-a To eliminate this limitation, three strategies were adopted in this study to characterize genes related to sepsis. First, recognized databases (HPO, TTD, TCMIP, DrugBank, DisGeNET, TCMSP) were used to search for sepsis correlated target genes (109 item genes).
Second, genes with significant transcriptional change in CLP group (Top 200 DEGs) were utilized to characterize targets related to sepsis.
Subsequently, WGCNA was conducted to identify co-expression modules and module-associated genes, defining sepsis-related targets (156). Thereafter, a multisource, comprehensive set of sepsis-related genes was obtained by integrating the aforementioned three sets of data, offering information fundamental for the analysis of network pharmacology. The "sepsis-related gene-MEL-putative target" interaction network was constructed using the well-supplied disease genes and drug putative targets. The hub genes scanned in this network have been partially identified (e.g., Saa3 and Fpr2) as potential targets of septic myocardial injury in previous studies. 31 N), and low-quality reads containing more than 20% of bases with qualities of <20. Hisat2 software was employed to align clean reads to mouse genome (version: mm10 NCBI). The SVA/Combat software 33 was utilized to alleviate any potential batch effects due to multiple sequencing runs using a raw count matrix ( Figure S1). The EdgeR package in R language was applied to normalize mRNA sequencing data. The DEGs between CLP and Sham groups were obtained using Deseq2 and the threshold was set as follows:

| GO terms and KEGG enrichment
To illustrate biological implications behind different gene sets including WGCNA core genes, DEGs, and hub genes, GO terms and KEGG pathway enrichment analyses were applied using the ClusterProfiler 34 package with cut-offs of p value < 0.05 in the Fisher's exact tests and an FDR-corrected p value < 0.1.

| Weighted gene co-expression network analysis
Quality-controlled count matrix was applied to identify significant gene modules and determining the core genes related to biometric, plasma and echocardiographic parameters of septic mice using the WGCNA package 35 in R. The filtering criteria for the expression data of the Sham and Sepsis group are as follows: (1) filter out the data whose sum of the expression levels of all samples is less than 10; (2) filter out the data whose coefficient of variation exceeds 50 in the sham group and sepsis group, respectively. The mouse biometric, plasma, and echocardiographic parameters were collected at 8 h post-CLP (Table S4). A weighted adjacency matrix was constructed based on the power adjacency a ij = jcor(x i ,x j )j β . 35 The cor operator denotes the Pearson's correlation coefficient in levels between gene i and gene j. The power of β (weighted coefficient) was set to be 8 (R 2 = 0.85) to establish a signed scale-free network. Next, a topological overlap matrix (TOM) was derived using the adjacency matrix data and the dynamic tree cut algorithm was performed to identify gene modules (which corresponded to colored branches of dendrogram). To decipher connections between gene modules and biometric, plasma and echocardiographic parameters of septic mice, Pearson's correlation coefficients between mouse traits and the module eigengene (ME, the first principal component of each gene module) were calculated. 35 In order to screen the core gene candidates, modules associated with the traits of interest, in our case "SV," "CO," "group," and "sepsis Subsequently, a total of 453 genes were collected as sepsis-related genes following integration and de-duplication of the three sets (Table S2).

| Network construction
The "sepsis-related gene-MEL-putative target" interaction network was analyzed using the public database STRING v11.5 24 between sepsis-related genes and putative targets of MEL. Official gene symbols of sepsis-related genes and MEL-putative targets were used as input to construct the network. The STRING parameter is 0.4 of confidence level and is 0.05 of FDR stringency. Network was visualized through Cytoscape software 42 (Version 3.5.1). To evaluate centrality, three topological characteristics 42 (degree, betweenness, and closeness) were calculated for each node in the network. For a given node, degree referred to the number of links to it, betweenness was the number of all shortest paths through it, and closeness corresponded to the inverse of the average shortest path length from it to all other nodes in the network. A node with larger values of the topological characteristics indicates a more important role in the interaction network. A hub gene was defined as a node whose degree, betweenness, and closeness values exceeded the median values of the overall nodes. KEGG pathway enriched by the hub genes was also displayed in the network.

| Rodent model of sepsis
The CLP model was used to induce sepsis as previously described 44 (Figure 4a,b). At 8 h post-CLP, the murine sepsis score (MSS) and anal temperature were measured. The MSS system involves observing seven components: appearance, level of consciousness, activity, response to stimulus, eyes, respiratory rate, and respiratory quality. Each of these variables is given a score between 0 and 4. 25 The established MSS score is the total of these seven components.

| Experimental design
MEL was dissolved in DMSO and given to mice every 2 days for 6 days before CLP injury. Mice were randomly assigned to the following groups: (a) mice injected with 1 ml/kg DMSO and subjected to sham surgery; (b) mice injected with 1 ml/kg DMSO and then subjected to CLP surgery; and (c) mice treated with MEL and then subjected to CLP surgery. In the third group, different levels of MEL (15, 30, or 60 mg/kg) were administrated to mice to monitor survival rate for 72 h, among which 30 mg/kg was selected for further analysis, including sepsis score, anal temperature, blood routine parameters, blood biochemical parameters, cardiac function, myocardial structure, and molecular signaling (Figure 4c).

| Echocardiography evaluation
Transthoracic echocardiography was performed using an animalspecific instrument (VisualSonics Vevo3100, VisualSonics, Toronto, ON, Canada) at 8 h post-CLP in mice. Mice were anesthetized with 3% isoflurane and 1 L/min 100% oxygen in an induction chamber for 1-2 min. Next, mice were, respectively, laid supine on a warm platform and kept anesthetized by 2% isoflurane until they lost body-righting reflex. Then, series of M-mode images F I G U R E 7 The mechanisms of MEL protecting against septic myocardial injury via inhibiting PIK3CG-related signaling pathway. MEL, melatonin at the level of papillary muscles were obtained. SV, CO, LVPWs, LVPWd, LVEDV, LVESV, and HR were measured using Vevo LAB 3.0.0. All measurements were based on three consecutive cardiac cycles.

| Histological staining
The myocardium was fixed in 4% paraformaldehyde and sectioned at a thickness of 4-5 μm. Morphological changes in myocardium were observed by hematoxylin-eosin (H&E) staining. The degree of myo-  Table S1. The reaction conditions were as follows: (1)  5.13 | Cell culture, treatment, and establishment of the LPS-induced cardiomyocyte injury model HL-1 cells were seeded and cultured as previously described. 45 The specific steps used to establish LPS-induced cardiomyocyte injury model were described below. First, HL-1 cells seeded in 60 mm culture dishes were allowed to grow to approximately 80% confluence.
Next, the medium was replaced with a fresh serum-free medium.
HL-1 cells were treated with MEL (3 h) before exposure to LPS (3 h) to examine the effects of MEL against LPS-induced cardiomyocyte injury. Finally, cell viability was detected using a Muse cell analyzer (Merck KGaA, Darmstadt, Germany). 5.14 | Cell viability and intracellular ROS generation assessments HL-1 cells were seeded and cultured as previously described. The cell viability and intracellular ROS generation were examined according to our previous study. 45

| Western blotting
Isolated cells were homogenized in RIPA buffer containing protease and phosphatase inhibitors (Beyotime Biotechnology, Shanghai, China).
Thirty-fifty μg of total protein extract was loaded to SDS-PAGE (8% or 10%) and transferred onto PVDF membranes (Millipore, Billerica, MA, USA), which were next blocked by 5% skim milk and incubated with pri-